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1  Introduction 


This  study  is  motivated  by  the  need  for  better  accuracy  in  atmospheric  turbulence  pre¬ 
dictions  for  the  upper  troposphere  and  lower  stratosphere  based  on  the  Penn  State/NCAR 
Mesoscale  Model  Generation  5  (MM5)  (Dudhia  1993;  Grell  et  al.  1994).  Although  tur¬ 
bulence  itself  is  not  captured  by  the  MM5,  MM5-predicted  atmospheric  profiles  of  winds 
and  temperatures  are  used  as  input  into  a  turbulence  parameterization,  such  as  the  Dewan 
model  (Dewan  et  al.  1993),  to  calculate  turbulence  along  paths  on  the  order  of  1000  km 
long,  providing  a  feasible  method  for  real-time  turbulence  predictions  for  use  in  aviation  and 
military  applications. 

Field  measurements  of  turbulence  and  other  atmospheric  quantities  were  made  to  assess 
the  accuracy  of  such  predictions  and  to  compare  them  to  turbulence  parameterizations  based 
on  observed  profiles  (Ruggiero  et  al.  2004).  This  revealed  cases  where  MM5-based  predictions 
failed  because  of  MM5  inaccuracies,  characterized  by  vertical  profiles  in  which  high-resolution 
features  were  lacking  or  severely  smoothed.  This  prompted  the  current  assessment  of  impacts 
on  MM5  accuracy  as  it  pertains  to  atmospheric  turbulence  calculations.  Gravity  waves  are 
of  interest  in  this  regard  because  gravity  waves  can  create  shear  instabilities  that  produce 
turbulence. 

Gravity  waves  are  the  focus  of  this  study  due  to  their  association  with  turbulence.  The 
MM5  captures  gravity  waves,  but  the  model  was  developed  as  a  weather-forecasting  tool  and 
its  simulations  are  subject  to  numerical  damping  that  was  not  optimized  for  gravity- wave 
forecasting.  In  the  context  of  turbulence  forecasting,  damping  of  gravity  waves  is  of  interest, 
especially  in  consideration  of  the  smoothed  MM5  profiles  in  cases  of  inaccurate  turbulence 
predictions. 

Related  research  by  Leutbecher  and  Volkert  (2000)  compared  several  MM5  simulations  to 
an  observed  gravity-wave  event  induced  by  flow  over  a  mountain.  This  showed  how  variations 
in  horizontal  grid  spacing  and  surface  friction  values  were  important  in  determining  the 
accuracy  of  the  simulated  wave  amplitude.  Horizontal  grid  spacing  affected  accuracy  because 
it  impacted  the  resolution  of  the  mountain  and  the  amount  of  artificial  dissipation,  which  is 
targeted  to  control  numerical  instability  of  waves  that  are  small  relative  to  grid  spacing.  A 
result  of  the  study  was  an  optimal  pairing  of  surface  friction  value  and  horizontal  resolution 
(the  finest)  for  the  given  wave  event. 

The  present  study  aims  to  separate  MM5  wave  propagation  by  the  model’s  dynamical 
core  from  other  model  functions  and  to  generalize  results  to  a  range  of  gravity  waves.  The 
dynamical  core  refers  to  that  part  of  the  model  that  is  adiabatic  and  invisid  and  thus  mainly 
pertains  to  advection  of  conserved  properties.  That  is,  propagation  is  modeled  without  the 
parameterizations  typically  used  to  represent  subgrid  scale  influences,  and  boundary  layer 
processes  so  that  the  effects  of  the  dynamical  core  are  isolated.  In  our  idealized  simulation, 
conditions  at  one  boundary  correspond  to  gravity  waves  entering  the  model  domain.  Then 
calculated  values  along  the  center  of  the  propagating  gravity  wave  are  used  to  evaluate  how 
the  amplitude  of  the  wave  evolves.  The  problem  simulated  has  an  analytical  solution  so  the 
observed  loss  in  wave  amplitude  can  be  assessed  in  terms  of  theoretical  results.  Further  to 
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generalize  the  simulation  results  we  conduct  a  numerical  analysis  of  the  linear  behavior  of 
the  discretized  model  equations. 

The  problem  simulated  is  the  two-dimensional  gravity-wave  problem  characterized  by  the 
X-shaped  wave  pattern,  known  as  St.  Andrew’s  cross,  shown  in  the  photographs  by  Mow¬ 
bray  and  Rarity  (1967).  A  density-stratified  fluid  with  constant  Brunt  Vaisalla  frequency, 
initially  at  rest,  is  disturbed  by  an  oscillator,  which  results  in  four  straight  gravity- wave 
beams  propagating  from  the  disturbance  location  and  forming  an  X  pattern.  The  nonlinear 
problem  is  solved  analytically  by  Tabaei  and  Akylas  (2003)  for  a  two-dimensional,  viscous, 
Boussinesq  fluid.  They  show  that  after  initial  transients  died  out,  their  solution  corresponds 
to  the  solution  of  Thomas  and  Stevenson  (1972)  who  validate  the  solution  with  experimental 
measurements.  It  turns  out  that  in  the  inviscid  case,  the  nonlinear  and  linear  solutions  are 
identical  and  prescribe  a  uniform  fluid  speed  along  each  beam  center.  Tabaei  and  Akylas 
showed  how  the  presence  of  viscosity  causes  the  beam-center  speed  to  decrease  as  the  beam 
propagates  away  from  the  disturbance. 

For  the  MM5  simulation  of  St.  Andrew’s  cross,  moisture,  rotation,  and  other  effects  that 
are  important  in  the  real  atmosphere  are  eliminated,  making  model  conditions  close  to  those 
of  the  Tabaei  and  Akylas  (2003)  solution.  The  compressibility  in  the  model  must  be  retained 
and  therefore,  the  theoretical  Boussinesq  results  for  the  beam-center  speed  are  adapted  to  a 
compressible  fluid  for  comparison  to  model  results.  Although  both  the  model  and  theory  are 
nonlinear,  small  disturbances  are  used  here  to  generate  a  linear  response  for  consistency  with 
numerical  analysis  described  later.  Results  from  eight  model  runs,  reflecting  different  grid 
sizes  and  wavelengths,  are  compared  to  the  theoretical  prediction  for  loss  in  beam  amplitude. 
This  shows  MM5  losses  can  be  greater  or  less  than  predictions,  depending  on  grid  and  wave 
parameters. 

A  linear  analysis  is  used  explore  the  sensitivity  of  the  MM5’s  numerical  dissipation  to 
changes  in  grid  and  wave  parameters.  It  is  noted  that  the  damping  in  the  MM5  is  strictly 
numerical,  as  real  viscosity  is  not  explicitly  simulated  by  the  dynamical  core  of  the  model. 
The  MM5  is  a  fully  compressible,  nonhydrostatic,  three-dimensional  model.  Variables  are 
updated  using  a  combination  of  explicit  time  stepping  with  artificial  dissipation  for  hori¬ 
zontal  momentum  equations  and  implicit  time  stepping  for  vertical  momentum  and  entropy 
equations.  The  terms  in  all  equations  are  partitioned  into  fast-changing  terms,  which  are 
updated  on  small  time  steps,  and  slowly-changing  terms,  which  are  updated  on  large  time 
steps. 

Numerical  experiments  result  in  observed  amplitude  changes  that  vary  with  grid  and  wave 
parameters.  The  observed  values  show  how  numerical  dissipation  in  the  MM5  dynamical 
core  compares  to  the  dissipation  that  would  be  present  with  molecular  viscosity.  Numerical 
dissipation  is  decreased  with  finer  wave  resolution  in  either  the  horizontal  or  vertical,  with 
greater  sensitivity  to  horizontal  resolution.  This  is  supported  by  both  numerical  experiments 
and  numerical  analysis.  Numerical  analysis  is  used  to  calculate  amplitudes  for  a  range  of 
horizontal  and  vertical  parameters,  showing  these  sensitivities. 
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2  Theory 


The  relevant  equations  and  their  analytical  solution  are  discussed  here  and  used  in  the 
next  section  where  the  numerical  method  is  analyzed.  The  two-dimensional  conservation 
equations  of  continuity,  momentum,  and  entropy  are  taken  for  a  compressible,  inviscid, 
adiabatic,  ideal  gas  and  linearized  about  a  motionless  background  state  with  constant  Brunt 
Vaisalla  frequency  and  temperature.  Although  approximate  solutions  adequately  describe 
gravity  waves  for  many  purposes,  the  compressible  case  is  most  relevant  since  compressibility 
is  an  essential  feature  of  the  MM5  numerical  method.  This  problem  leads  to  a  gravity-wave 
solution  and  the  associated  dispersion  relation.  When  the  problem  is  further  simplified 
using  the  Boussinesq  approximation,  a  gravity-wave  beam  solution  is  obtained.  It  is  shown, 
by  Tabaei  and  Akylas  that  the  Boussinesq  solution  for  the  linearized  problem  is  the  same  as 
that  for  the  nonlinear  problem.  Furthermore,  they  give  the  solution  for  the  Boussinesq  case 
with  viscosity.  All  of  these  solutions  are  considered  here. 

The  perfect  gas  law, 

P  =  pRT  (1) 

where  p ,  p,  and  T  are  pressure,  density,  and  temperature,  respectively,  and  where  R  is  the 
gas  constant  for  air,  is  linearized  using  the  sum  of  a  background  state  quantity,  indicated  by 
the  subscript  0,  and  a  perturbation  quantity, 

p  —  Pq -jr  p  +  ■  ■  •  ,  p  =  po  +  p  +  •  •  •  and  T  =  To  +  T  +  •  •  •  .  (2) 

Requiring  that  Equation  (1)  holds  for  the  background  quantities  and  ignoring  terms  of  second 
order  leads  to  the  linearized  perfect  gas  equations, 


P  —  PqTR  +  pToR. 


(3) 


Since  it  is  already  linear  the  hydrostatic  relationship  for  the  perturbation  quantities 


dp 

Tz 


=  ~P9 


(4) 


has  the  same  form  as  for  the  full  and  background  variables. 

The  background  state  sound  speed,  c0,  and  Brunt  Vaisalla  frequency,  A,  are  constants 
in  this  treatment,  given  by 


Co  =  7RT0,  and  N  =  — a/t  -  R  (5) 

Co 

where  7  is  the  ratio  of  the  specific  heat  of  dry  air  at  constant  pressure  to  the  specific  heat 
of  dry  air  at  constant  volume  and  g  is  the  gravitational  constant. 

In  preparation  for  the  numerical  analysis,  the  governing  equations  are  shown  in  the  form 
used  in  the  MM5  and  without  the  added  damping  terms  from  the  MM5.  The  momentum, 
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pressure,  and  temperature  equations  are 


du  1  dp 
dt  po  dx 

0, 

(6) 

dw  1  dp 

p 

(7) 

dt  ^  po  dz 

- 9 — , 

Po 

P0Cq<5  -  p0gw  = 

0, 

(8) 

dT 

1  dp  g 

(9) 

-  = 

- - - 

dt 

Po  &p  dt  Cp 

5  = 

du  dw 
dx^  dz' 

(10) 

Equations  (6)  -  (10)  match  those  shown  in  Grell  et  al.  (1994,  referred  to  hereafter  as 
G94),  when  the  MM5  mapping  factor  is  m  =  1,  the  reference  pressure  p*  is  constant  (as  it 
would  be  for  flat  terrain),  and  using  2  as  a  vertical  coordinate.  Specifically,  Equations  (6)  - 
(9)  above  can  be  obtained  respectively  from  Equations  (2.2.1)  or  (2. 5. 1.1),  Equations  (2.2.3) 
or  (2. 5. 1.3),  Equation  (2.2.4),  and  Equation  (2.2.5)  in  G94. 

Equations  (6)  -  (10)  are  equivalent  to  those  presented  in  Lighthill  (1978)  for  the  flow  of  an 
inviscid,  linearized,  compressible,  constant  N  fluid.  The  momentum  equations  above,  Equa¬ 
tions  (6)  and  (7),  are  both  represented  by  the  (vector)  Equation  (16)  in  Lighthilbs  Section 
4.1.  The  pressure  equation,  Equation  (8)  above,  can  be  obtained  by  combining  Lighthill’s 
equation  for  a  reversible  process,  or  Equation  (33)  in  Section  4.1,  with  the  continuity  equa¬ 
tion,  or  Equation  (30)  from  Section  4.1.  The  temperature  equation,  Equation  (9)  above,  can 
be  derived  from  Equation  (33)  in  Lighthill  Section  4.1,  along  with  the  linearized  perfect  gas 
law,  the  hydrostatic  relation,  and  the  relations  between  constants  in  Equation  (5).  Thus, 
the  solution  from  Lighthill  applies  here. 

The  solution  to  Equations  (6)  -  (10)  takes  the  form 


p  =  TZ  exp  (icj)),  p0u  —  Uexp{i(p),  p0w  =  Wexp(i^),  p  =  V  exp(i(f>),  (11) 


where 


4>  =  kxx  +  kzz  —  ut , 


(12) 


uj  is  the  temporal  frequency,  and  kx  and  kz  are  the  spatial  frequencies  in  the  x-  and  2- 
dimensions,  respectively.  This  leads  to  the  dispersion  relation 


to 

co 


kt  +  k2  —  ikz 


7  N2 

(7  -  1)3. 


u 


+  kiN2  =  0, 


(13) 


as  shown  in  Lighthill  (1978).  Provided  that  the  background  density  varies  slowly  over  the 
length  of  a  gravity  wave,  the  approximate  dispersion  relation  is 


u> 


2 


N2k2x 

%  +  k? 


(14) 
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which  coincides  with  the  Boussinesq  dispersion  relation.  The  wave-dependent  coefficients  in 
(11)  are  related  by 

V  uj  V  N2  —  u2  TZ  iui  ikz  V  .  . 

«  =  v  w  =  w  =  (  > 

The  nonlinear,  two-dimensional,  Boussinesq  beam  problem  was  solved  by  Tabaei  and 
Akylas  (2003)  for  a  stably  stratified  fluid  with  background  density,  po,  that  varies  exponen¬ 
tially  with  height.  The  gravity-wave  dispersion  relation,  which  coincides  with  the  linear  case, 
is 

(lj/N)2  =  sin2  9,  (16) 

where  6 ,  the  wave  propagation  angle  in  the  x-z  plane,  is  defined  by 

(kx,kz)  =  K  (sin  8,  —  cos  6) .  (17) 


Wave  groups  move  with  the  group  velocity  given  by 


^  N  cos  8 1  . 

C  g  =  — — — (cos£/,sm0). 


(18) 


Provided 

K  »  (19) 

Po 

these  results  are  a  good  approximation  to  the  linear,  compressible  solution,  which  is  presented 
later  and  can  be  found  in  Lighthill  (1978). 

Along  the  center  of  a  beam,  the  fluid  speed  is  constant  for  the  inviscid,  Boussinesq  case. 
Tabaei  and  Akylas  (2003)  showed  that  when  viscosity  is  accounted  for,  beam-center  speed, 
V,  away  from  the  disturbance  decays  according  to 


V  oc 


1 

IF5’ 


(20) 


where  £  is  the  beam-following  coordinate, 


£  =  x  cos  9  +  zsinff. 


(21) 


The  beam  amplitude  prediction  is  reinterpreted  here  for  a  compressible  fluid,  for  consistency 
with  MM5  assumptions.  In  a  compressible  fluid,  the  factor  of  y£oo  constitutes  the  leading 
correction  to  Boussinesq  results  (Lighthill  1978).  Therefore,  the  beam  amplitude  in  the 
compressible,  viscous  case  used  here  is 


Vy/po  oc 


1 

IF5' 


(22) 


This  amplitude  prediction  is  used  to  gauge  whether  MM5  beam  amplitudes  are  diminished 
by  MM5  damping  more  or  less  than  they  would  be  by  molecular  viscosity. 
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3  Numerical  Analysis 


To  examine  how  damping  of  gravity  waves  in  MM5  varies  with  grid  and  wave  parameters, 
the  discrete  forms  of  (6)  -  (9)  are  analyzed  using  the  MM5  numerical  scheme.  Many  of  the 
techniques  used  here  are  shown  in  Durran  (1999).  A  result  is  a  numerical  dispersion  relation 
from  which  the  wave  frequency  and  the  wave  amplification  factor  are  found. 

The  MM5  numerical  scheme  is  based  on  a  form  of  Equations  (6)  -  (10)  with  added 
damping  and  averaging  to  accommodate  the  staggered  grid.  The  equations  are  rewritten 
with  these  features  and  using  .ft  =  cp  -  <q,,  7  =  cp/ct,,  and  c2a  =  'ypo/po .  The  result,  after 
some  rearrangment,  is 


8u  8  f  p  \  86 

dt  0  dx  ^  poco )  1  dx 


dw 

dt 


+  Cq 


8  (  p 


8z  \poco 
d  f  V 
dt  VPoCo 

JL(2sL\ 

dt  V  Co  ) 


Co 


P0C0 


+  cq6  -  —w 
Co 


d_ 

dt 


P 

P0C0 


D(u), 

g{ 7  - 1)  f  cPT  p 


Co 


Co  P0C0 


0, 


°-W-dU- 

co  \PoCq 


+  D 


+  D(w), 


(¥) 


(23) 

(24) 

(25) 

(26) 


where  the  D  terms  are  fourth  order  diffusion  terms,  the  a  term  is  divergence  damping,  and 
the  overbar  indicates  vertical  averaging.  The  damping  terms  are  used  in  MM5  for  numerical 
stability  and  do  not  reflect  a  physical  process. 

The  vertical  coordinate,  2,  used  here  differs  from  the  pressure-based  coordinate  used  in 
the  MM5,  since  use  of  the  MM5  coordinate  would  unnecessarily  complicate  the  analysis.  For 
the  given  base  state  the  results  would  not  change  since  pressure  levels  and  2  levels  coincide. 
In  the  MM5  simulations  the  pressure  levels  are  chosen  to  correspond  to  equal  increments  in 
2. 

The  MM5  equations,  approximated  by  the  linearized  equations  in  (23)  -  (26)  with  (10), 
are  discretized  with  centered  differences  in  space.  Note  in  what  follows  that  the  MM5 
staggers  the  grids  of  the  different  variables.  In  the  horizontal  domain  the  Arakawa  D  grid 
is  used  in  which  the  thermodynamical  variables  (temperature  and  moisture)  and  vertical 
velocity  are  defined  at  the  full  grid  points  (with  integer  indices)  and  the  horizontal  wind 
components  are  defined  at  the  half  grid  points,  he.,  at  the  corners  of  the  grid  boxes  centered 
on  the  full  grid  points.  In  the  vertical,  vertical  velocity  is  defined  at  the  full  grid  points  and 
all  other  variables  are  defined  at  the  half  grid  points.  Therefore,  the  first  derivatives  of  a 
variable,  v,  accounting  for  grid  staggering,  are  approximated  in  the  MM5  by 


dv  ^  Vj+i/2  —  Vj- 1/2  A  8v  Wfc  +  1/2  —  Vk- 1/2 

— —  ^  and  —  ~ 

OX  Xj+l/2  ~  Xj-l/2  <V2  Zk+ i/2  —  Zk-l/2 


(27) 


and  the  vertical  average  is 


v  — 


Vk  +  l/2  +  Vk- 1/2 
2 


(28) 
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For  a  wave  mode  represented  by 


v  =  V  exp  (ikxx  +  ikzz)  exp  (29) 

and  using  Ax  =  Xj+ 1/2  —  ^-1/2  and  Az  =  2^1/2  —  z*,- 1/2,  the  approximate  derivatives  and 
the  vertical  average  can  be  written, 


dv  _  isin  kxAx/2 

8z~d‘V=  Ax/2 

(30) 

3u  2  sin  kzAz/2 

az K  =  Az/2 

(31) 

U  =  azv  =  coskzAz/2v. 

(32) 

An  artificial  fourth  order  diffusion  is  included  in  the  MM5  to  selectively  dissipate  the 
smallest  scale  waves.  Otherwise,  because  of  aliasing,  significant  energy  would  accumulate 
in  these  scales.  The  fourth  order  diffusion  terms  are  defined  so  that  D(v)  approximates 
D4:div/dxi,  with  D4  a  constant  dependent  on  grid  parameters.  This  approximation  is  ac¬ 
complished  by 


D(v) 

Da 


D 4 


Vj+ 2  -  +  6vj  -  4vj-i  -  Vj- 2 

Az4 


.003 


Az4 
At  ’ 


(33) 

(34) 


where  the  Z?4  value  is  approximated  using  Section  5.1  in  G94,  recognizing  that  the  MM 5 
coefficient  called  Ku  is  dominated  by  the  Kjj 0  contribution  for  gravity  waves  and  At  is  the 
so-called  small  time  step.  Using  the  wave  mode  given  above,  the  fourth  order  diffusion  term 
is 

012 

D{v)  =  Dxv  =  —  — —  (cos  kx Ax  -  l)2  v.  (35) 

At 

Discretizing  the  spatial  derivatives  and  using  the  above  notation  for  averages,  Equations 
(23)  -  (26)  and  (10)  are  written  in  differential-difference  form  as 


du 

dt 


+  cq9x 


P 

Poc 0 


&xdxS 


—  Dxu, 


dw 

dt 


VP0C0/  cq  \P0C0J  cq  \  cq  PoCo/ 


,Poc0 

s!_(  _p 

dt 


d_ 

dt 


PoCo 

CpT 

Co 


+  CoS  —  — w 
Co 


=  0, 


d_ 

dt 


V 

PoCo 


9-  n 
— w  —  D 
Co 


+  D  SeL 

x  1  *-/a:  ) 

PoCo  Co 


where  <5  is  now  interpreted  in  terms  of  the  difference  notation  so 


5  =  dxu  +  dzw. 


(36) 

(37) 

(38) 

(39) 

(40) 
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The  above  equations  are  next  subject  to  the  time  differencing  for  the  time  derivatives, 
and  time-level  assignments  for  all  variables.  In  summary,  the  MM5  scheme  partitions  terms 
in  the  momentum  and  pressure  equations  into  fast-changing,  or  acoustic,  terms  and  slowly 
changing  terms.  The  fast-changing  terms  are  updated  on  small  time  steps,  At,  while  the 
slowly-changing  terms  are  updated  on  large  time  steps,  2 At.  In  addition,  some  of  the  slowly- 
changing  terms  are  updated  at  times  t  -  At,  t  +  At,  etc.,  while  others  are  updated  at  times, 
t,  t  +  2 At.  So  at  any  given  time  step,  multiple  time  levels  must  be  represented.  The  scheme 
is  described  in  more  detail  below  and  used  to  rewrite  the  differential-difference  equations  as 
difference-only  equations.  This  combination  of  split-explicit  and  odd/even  centered  (or  “leap 
frog”)  time  steps,  along  with  a  time  filter  to  keep  the  odd  and  even  solutions  from  diverging, 
is  used  in  MM5  to  simultaneously  optimize  the  efficiency  and  accuracy  of  the  model. 

The  fast-changing  terms  in  the  MM5  appear  on  the  left  hand  side  of  (36)  -  (38)  and  slowly- 
changing  terms  appear  on  the  right  hand  side.  The  slowly-changing  terms  in  the  vertical 
momentum  equation  (37)  are  individually  susceptible  to  acoustic  oscillations,  but  their  sum 
is  equivalent  to  the  slowly-varying  potential  temperature  term,  p0gQ/Qo,  where  0o  and  0  are 
the  leading  terms  in  the  expanded  potential  temperature  and  ©o  =  Po/(po7?)(ps/po)^7-1/7\ 
with  a  reference  pressure  ps . 

The  acoustic  terms,  which  appear  on  the  left  hand  side  of  (36)  -  (38),  are  updated  at 
every  small  time  step,  At,  while  the  other  terms  are  updated  on  large,  “leap  frog”  (i.e., 
centered)  time  steps  of  size  2 At,  with  4Ar  =  2 At,  typically.  The  small  time  steps  involve 
stepping  forward  in  time,  so  that  for  u  for  example, 

du  _  u{t  +  At)  —  u(t)  . 

8t~  At  ’  ^ 

and  likewise  for  w  and  p/(poCo).  In  MM5  the  u  equation  is  stepped  forward  explicitly  and 
then  the  updated  u  values  are  used  in  the  w  and  p/(poco)  equations.  Equations  (37)  and 
(38)  are  made  implicit  in  the  small  time  steps  by  using  time-averaged  values  for  pressure 
and  vertical  velocity.  As  described  in  Section  2.5.1  in  G94,  occurences  of  p/(poCo )  among 
the  fast-changing  terms  in  the  w  equation,  or  on  the  right  hand  side  of  Equation  (37),  are 
replaced  by  the  time-averaged  value 


1  +  P  V 

2  PqCq 


+  At)  + 


P  ,  s 
2  PoQ) 


(42) 


where  0.2  <  j3  <  0.4.  Similarly,  occurences  of  w  on  the  right  hand  side  of  Equation  (38)  are 
replaced  by  the  analogous  time  averaged  value  of  w.  For  analysis  purposes  the  median  value 
P  =  0.3  is  adopted. 

The  temperature  is  updated  on  leap  frog  steps,  which  leap  from  t  —  At  to  t  +  At.  The 
nonacoustic  terms,  other  than  the  fourth  order  diffusion  terms,  are  evaluated  at  the  large 
time  levels  in  between  these  leaps,  or  t,  t  +  2 At,  etc.  The  fourth  order  diffusion  terms  are 
evaluated  at  t  —  At,  t  +  Af,  etc.  To  write  the  difference  equations,  time-level  superscripts 
are  used,  with  m  and  n  superscripts  representing  small  and  large  time  steps,  respectively. 
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The  small-step  difference  in  (41)  is  written 


u(t  +  At)  —  u(t)  um+1  —  um 


At 


At 


(43) 


while  the  large-step  difference  to  approximate  the  time  derivatives  in  Equation  (39)  use  the 
form. 


££^pn+l  £z'jpn-l 

co _ co _ 

At 


(44) 


Starting  from  the  step  n  -  1  =  m,  the  first  small  time  step  gives  the  difference  equations 


u 


m+ 1  _  m 


W 


P 


m+ 1 


PoCo 


Co 


u 
+Ar 


-cq8x  ( )  +  axdx(dxum  +  dzwm)  +  Dxun  1 

\poCo  1 


m+ 1  _ 


-  At  La,  +  ' ) 

V  Co  J  \  2  p0co  2  pqCqJ 

+A  Tazdz(dxum+1  +  8zwm) 

+AT£(mliI!,  {<£1  _  jq + 

Co  \  Co  PoCo/ 

—  -  Arcoaxum+1 
PoCo 

+Ar  f-Co52  +  —  az^  ( +  ^Awm 


CpTn+1  CpT 


co 

in-l  »!«+! 


CO 


+ 


P 


PoCo  PoCo 


,n— 1 


CpT' 


-At  (  -  Dx— — -  +  Dx 

Co  PoCo  Co 


in— 1 


(45) 


(46) 


(47) 


(48) 


Note  that  with  four  small  steps  for  every  leap  frog  step,  every  other  m-based  step  coincides 
with  an  n-based  step,  so  that  m  =  n  —  l,m  +  2  =  n,  and  m  +  4  =  n  +  1.  The  distinction 
between  m  steps  and  n  steps  in  Equations  (45)  -  (48)  is  retained  because  equal  values,  such 
as  m  and  n  —  1,  for  the  step  represented  above,  will  not  be  equal  in  the  subsequent  small 
step  when  m  is  updated  to  m  + 1  and  n  —  1  is  held  constant.  Note  Equations  (45)  -  (48)  are 
valid  for  any  number  of  small  time  steps  per  big  time  step.  These  equations  will  be  used  to 
develop  a  set  of  equations  for  the  specific  case  of  four  small  time  steps  for  a  leap  frog  step. 

Equations  (46)  and  (47)  are  rearranged  to  eliminate  occurences  of  m  +  1  level  quantities 
on  the  right  hand  side.  Specifically,  Equations  (46)  and  (47)  are  solved  simultaneously  to 
eliminate  wm+1  and  pm+1/(p0c0)  and  Equation  (45)  is  used  to  replace  occurences  of  um+1. 
Rewriting  the  u  equation  along  with  the  result  of  the  rearranged  w  and  p/(poCo)  equations 
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gives 


Here 


,m+ 1 


W 


,m+l  _ 


P 


.m+1 


PoCo 


+  A™wm  +  P 

+  Anu-lun- 

-1 

) 

(49) 

Poco 

£>m  +  £(>m  +  B"1  pm 

PU 

+  Hpn— 

pn 

+  B?^— 

(50) 

Poco 

iBrv-'  +  srv-1, 

pPoCo 

PoO) 

c>m  +  c;>m  +  c?  pm 

vn 

+  c$— 

„n 

(51) 

p  Pa  Co 

+  Cn~lun~l  +  Cn-lwn-\ 

v  Pqc0 

PoO) 

A™ 

A™ 

a: 


1  +  axdx, 

(5X  dx  dz , 

— ^  4n~1  —  n 


(52) 

(53) 

(54) 


j^m+l 


Bl 


1 +  («;  -  ®  i'i±^V 


R?  = 


~2  2  )  ’ 

{az  +  4)  (^ir)  (1  +  §jdx)dx 

Rm+1  ’ 

1  -  (5x  +  dz)  (^)  f (a,  -  d2)  (i=£)  +  did. 


a-r 


Qm+1 


B™  =  - 


(a,  +  a)  (i  +  ^4) 


£Jm+l 


£? 


Bp  = 


b: 


n— 1 


(7  -  IK 

Dm+1  5 

(7  -  l)flz 

Rm+1  5 
w 

(az  +  dz){^)dxDx 
5m+1 


(55) 

(56) 

(57) 

(58) 

(59) 

(60) 
(61) 
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and 


C? 

=  -4  +  (a,  -  4)  (li^  w  - 

(62) 

cr 

=  (a,  -  4)  -  4&S., 

(63) 

c. m 

=  1  +  d2x  +  (az-dz)^-B™, 

(64) 

'—/rp 

=  (az-dz)^BZ, 

(65) 

c; 

=  (a,-4)~B;, 

(66) 

cr1 

(67) 

/nm— 1 

=  (5,-4)ii£Br', 

(68) 

with  the  dimensionless  quantities, 

az  =  ^—~az,  dx  =  c0Ardx,  dz  =  c0Ardz,  and  ax  -  (69) 

Cq  At  Cq 


defined  for  convenience. 

The  discrete  system  is  appended  with  the  following  equation, 

CpTn  _  CpTn+1  OpTn~l 

c0  2c0  2c0  ’ 


(70) 


which  simply  averages  temperature  values.  Note  that  all  the  coefficients  needed  to  describe 
the  linearized  MM5  equations  are  ultimately  only  dependent  on  the  wave  numbers  kx  and 
kz]  the  “delta’s”  Ax,  A z,  At,  and  At;  the  physical  constants  g,  Cp,  and  7;  the  basic  state 
parameters  Co  and  po\  and  the  pressure  term  averaging  parameter  (3. 

Equations  (45)  -  (48)  along  with  Equations  (49)  -  (70)  provide  the  rules  for  discrete 
time  stepping  for  any  number  of  small  time  steps  per  large  time  step.  To  continue  the 
analysis  we  must  now  specify  how  many  small  time  steps  correspond  to  a  leap  frog  time 
step.  The  standard  MM5  always  defaults  to  four  small  time  steps  per  leap  frog  step,  and 
this  arrangement  is  used  in  our  simulation  experiments.  Thus,  Equations  (49)  -  (51)  are 
written  for  m  +  2,  m  +  3,  and  m  +  4,  by  advancing  the  m-based  variables,  while  leaving 
the  n-based  variables  at  the  levels  specified.  For  example,  using  Equation  (50)  to  write  the 
expression  for  wm+2  gives 

wm+ 2  =  B™+lum  +  B™wm+1  +  B™- - +  B£^—  +  BZ-Z—  (71) 

p  P0C0  p  P0C0  1  P0C0 

+B£- V-1  +  BZ-'wn-\ 
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Note  that  the  coefficients  B ™  does  not  actually  depend  on  u  or  m.  The  super-  and  sub¬ 
scripts  were  defined  to  match  the  variable  um  that  B ™  multiplies.  This  is  true  for  all  the  A , 
B,  and  C  coefficients. 

Writing  Equations  (49)  -  (51)  for  the  four  small  time  steps  results  in  twelve  equations, 
which  are  combined  to  reduce  the  number  of  equations  for  u,  w,  and  p/{poCo)  to  six.  Before 
combining  equations,  the  levels  m,  m  +  2,  and  m  +  4  are  renamed  using  the  coinciding  n- 
based  levels,  n  —  1,  n,  and  n  +  1,  respectively.  The  remaining  m-based  levels,  m  +  1  and 
m  +  3,  are  eliminated  using  substitution.  For  example,  in  the  equation  for  um+2,  now  called 
un,  occurences  of  um+1  and  pm+1/(Poco)  are  replaced  using  the  equations  for  those  variables. 

The  six  equations  for  u,  w,  and  p/(pqCq)  are  combined  with  the  two  equations  for  cpT/cq , 
namely  Equations  (48)  and  (70),  for  a  system  of  eight  linear  equations  with  eight  unknowns, 
containing  stricly  n-based  time  levels.  To  write  these  equations  compactly,  let 

i/>n=(V  wn  (72) 

^  \  poco  co  J  '  v  > 

so  the  eight  equations  are 

r+l  =  Lur  +  Lur-\  (73) 

4>n  =  L21iB  +  W71-1.  (74) 


0 

0 

0 

0 

B™B»+B™C2+B% 

0 

0 

wsM^ssimsm 

•5  (C?A%  + 

C%B™  +  C™C™) 

•5  (C£M£  + 

c;b:  +  c;c™- 

— a  ) 

_ 

■5  (C™A™  + 

C™(B™  +  Bp  + 

.5  (C"BJ+C”C?+ 

c?) 

and  the  elements  of  L2 2  are 


Anu~\A™  + 1)  + 
(A™)2  +  A™(B™  + 
Bl~l)  +  A™(C™  + 

cr1) 

A%(A2  +  B%  + 
BZ-')  +  A™(CZ  + 

cr1) 

0 

BZ{A™  +  Anu~l)  + 
B™(B™  +  B"_1)  + 
S-(C-  +  Crx)  + 

sr1 

B™A™  +  fir1  + 
B™(B™  +  5rx)  + 

£Pm(c™  +  qr1) 

DmAn 

0 

Now  the  numerical  scheme  is  reduced  into  the  steps  in  Equations  (73)  and  (74),  which 
can  be  thought  of  as  a  sequence  of  odd  and  even  steps.  To  distinguish  these  steps,  the  levels 
are  renamed  so  that 


iP2n  =  Lai^  +  W71-1,  (75) 

■02n+1  =  iui^  +  W71"1-  (76) 

These  equations  completely  describe  the  numerical  scheme,  excluding  the  time  filter. 

With  the  Asselin  (1972)  time  filter,  the  leap  frog  step  in  Equation  (76)  is  changed  to 

V,2n+1  =  In  +  Ln^ 7,  (77) 

where  the  time  filter,  indicated  by  the  overbar,  is  given  by 

ip2n  =  (1  —  2v)ip2n  +  i/i/)2n+1  +  vip2n-1,  (78) 

with  v  =  0.1  as  the  MM5  filter  constant.  Combining  Equations  (75)  and  (77)  gives 

iP2n+1  =  LnAil>2n-1  +  Ln ^7,  (79) 


where,  in  terms  of  the  identity  matrix  J, 


A  =  (/-L21)-1L22,  (80) 

so  the  scheme  is  completely  described  by  (79)  with  the  filter  defined  by  (78). 

In  order  to  obtain  a  single  equation  for  the  scheme,  in  which  filtered  variables  are  rep¬ 
resented  in  terms  of  unfiltered  variables,  the  filter  is  applied  to  both  sides  of  Equation  (75), 
which  when  solved  gives 

ijj2n  =  A'i/>2n_1,  (81) 
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using  the  property  that  Aip  —  Apj.  The  filter  in  Equation  (78)  can  now  be  rewritten  as 

Aip2n~l  =  (1  -  2u)A,ip2n~1  +  vtfn+1  +  (82) 

Solving  for  i/’2"-1,  this  gives 

=(A-  Iv )_1  [(1  -  2i')Ai/>2n_1  +  i^2n+1]  .  (83) 

This  is  used  to  eliminate  the  filtered  variable  in  Equation  (79)  giving  the  final  linear  evolution 
equation 

ip2n+1  =  Aip271-1  (84) 

where 

A=[I-  vLn{A  -  viy1] [LnA  +  (1  -  2 v)L12{A  -  uiylA]  .  (85) 

Next,  it  is  shown  how  this  result  is  used  to  obtain  the  amplitude  factor  and  the  numerical 
frequency.  Assuming  solutions  of  the  form 

'j/j2n+1  =  ip0  exp  {ikxx  +  ikzz)  exp  (— iujt2n+1),  (86) 

where  ipo  is  the  constant  vector  (uq  w0  Po/{poCq )  cpTq/co)t  that  we  must  determine. 
Equation  (84)  can  be  written 

0  =  [.A  —  1  exp  (— iu>2At)}  tpo  exp  {ikxx  +  ikzz)  exp  (—iu(t2n  —  At)),  (87) 

0  =  [A  -  IX]  (88) 

using  A  =  exp  (—iu2At).  Nontrivial  solutions  for  'ip2n~l  exist  for 

det  (A  -  IX)  =  0,  (89) 

which  leads  to  four  eigenvalues  (values  for  A),  from  which  the  appropriate  value  for  the 
gravity  wave  of  interest  is  selected,  as  discussed  below.  Each  eigenmode  then  simply  evolves 
according  to 

ip2n+1  =  Xip2n~l ,  (90) 

that  describes  the  effect  of  the  numerical  scheme.  The  amplitude  and  frequency  information 
are  contained  in  the  real  and  imaginary  parts  of  A,  or  in  terms  of  frequency,  they  come  from 


A  =  exp  (— i(ujT  -f  iuji)2At)  =  exp  (-iur2At)  exp  (ui&At).  (91) 

So  over  a  leap  frog  step,  the  amplitude  of  the  gravity  wave  of  interest  is  adjusted  by  a  factor 
of 


|  A  |  =  exp  (u>i2At)  =  Re(A) 
and  the  numerical  frequency  of  the  wave  is 


(92) 


~  Re(  2 At  )' 


(93) 
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D 

c. 

(+>+) 

+ 

(+>+) 

(+>■) 

(+>+) 

(■>+) 

■SSI 

(+.+) 

(+r) 

(">+) 

.  + 

(">+) 

(+>+) 

E 

(-i+) 

B 

(+)+) 

(+,-) 

Table  1:  Association  between  frequency  sign  and  gravity  wave  travel  direction.  Here  C  —  g 
is  the  group  velocity  and  Cp  is  the  phase  speed. 


Solving  Equation  (89)  leads  to  four  eigenvalues,  from  which  the  one  associated  with  an 
upward-traveling  gravity  wave  is  selected.  The  selection  is  based  on  the  magnitude  and  sign 
of  the  frequency.  Two  eigenvalues  lead  to  high  frequencies,  on  the  order  of  0.08  s-1,  which 
imply  oscillation  periods  on  the  order  of  1.3  minutes.  These  frequencies  are  higher  than 
the  Brunt  Vaisala  frequency  of  0.012s-1  and  are  therefore  discarded  as  describing  gravity 
waves.  The  other  two  eigenvalues  lead  to  frequencies  on  the  order  of  0.0008  s-1,  which  imply 
oscillation  periods  on  the  order  of  140  minutes,  which  are  consistent  with  gravity  waves. 
This  leaves  two  gravity  wave  eigenvalues,  which  are  opposite  in  sign  and  nearly  equal  in 
magnitude.  These  are  the  two  beams  of  the  cross. 

The  dispersion  relation  and  group  speed  from  the  idealized  case,  given  in  Equations  (16) 
-  (18),  can  be  used  as  a  guide  to  associating  frequency  sign  with  gravity  wave  direction 
of  travel.  For  example,  if  the  wavenumber  (kx,  kz)  signs  are  (+,-)  and  the  frequency  (to) 
is  positive,  then  the  group  speed  vector  ( Cg )  has  signs  (+,+)•  The  other  combinations  are 
given  in  Table  1  for  use  in  interpreting  calculations  in  association  with  the  numerical  analysis 
as  described  in  the  following  sections. 

Numerical  dissipation  is  typically  quantified  in  terms  of  the  order  r  in  A  =  1  —  Cfr , 
where  C  is  a  constant,  and  /  depends  on  grid  and  wave  parameters  (Strikwerda  1989).  In 
the  current  problem  /  is  not  readily  determined  and  therefore  the  quantity  V,  which  varies 
like  r\nf,  is  used.  The  amplification  factor  is  related  to  V  by 

In  (1  —  Aa)  =Va  +  Ba,  (94) 

where  Aa  =  Ay  ,  Ba  is  a  constant,  and  the  subscript  a  refers  to  results  based  on  analysis. 
The  constant  Ba  is  combined  with  its  counterpart  for  experimental  dissipation  as  discussed 
in  Section  5. 
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Figure  1:  MM5  domain  and  boundary  disturbance. 


4  Numerical  Experiment 

The  MM5  is  implemented  for  an  idealized  atmosphere  for  eight  sets  of  grid  and  wave 
parameters.  Simulations  are  in  a  single  domain  without  moisture,  Coriolis  effects,  a  plan¬ 
etary  boundary  layer,  or  tropopause.  The  terrain  is  uniform  and  the  background  state  of 
the  atmosphere  is  given  by  an  ideal  gas  in  hydrostatic  balance  with  no  flow  and  a  tempera¬ 
ture  profile  of  the  standard  atmosphere  (Holton  1992),  resulting  in  a  nearly  uniform  Brunt 
Vaisala  frequency  of  0.012s-1.  Such  conditions  should  be  steady,  but  because  of  boundary 
conditions  and  the  effects  of  descretization  the  specified  state  is  not  exactly  steady.  There¬ 
fore  a  preliminary  forecast  is  run  with  constant  boundary  conditions  to  allow  the  model  to 
make  any  small  adjustments  required  to  reach  a  steady  state. 

A  two-dimensional  disturbance,  satisfying  the  gravity- wave  solution  in  (11)  -  (15)  is 
introduced  at  one  lateral  boundary,  approximately  centered  in  the  model  vertical  domain. 
The  disturbance  has  a  vertical  extent  equal  to  one  vertical  wavelength,  above  and  below 
which  the  disturbance  decays  as  shown  in  Figure  1.  The  other  three  lateral  boundaries  have 
no-flow  conditions.  In  order  to  simulate  the  two-dimensional  problem  in  a  three-dimensional 
atmosphere,  the  disturbance  is  applied  uniformly  across  the  third  dimension  (the  N-S  or  y 
dimension  in  MM5  notation). 

The  resulting  flow  contains  an  upward-traveling  beam  and  a  downward-traveling  beam, 
forming  the  right  half  of  St.  Andrew’s  Cross.  Once  the  beams  have  developed,  cross  sections 
involving  the  N-S  dimension  are  examined  to  confirm  qualitatively  that  the  flow  is  uniform 
in  the  N-S  dimension,  indicating  two-dimensional  flow.  Cross  sections  showing  half  of  the 
X-pattern  are  examined  at  several  times  to  confirm  that  wave  quantities,  including  wave 
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Figure  2:  Example  of  MM5  gravity-wave  response  as  shown  in  a  cross-section  of  the  pertur¬ 
bation  density  field. 


numbers,  frequencies,  beam  angle,  and  group  and  phase  velocities,  are  consistent  with  theory. 

Reflections  from  the  top  of  the  domain  are  observed  and  found  to  interfere  eventually 
with  the  waves  of  interest.  To  exclude  the  influence  of  these  reflections,  a  final  time  is 
selected  for  each  wavenumber  so  there  is  sufficient  time  for  beams  to  form,  yet  not  enough 
time  for  reflections  to  have  a  noticeable  influence.  Beam-center  velocities  in  upper  beams 
are  measured  from  MM5  data  and  combined  with  the  known  background  density  to  obtain 
beam-center  amplitudes  used  in  dissipation  evaluations. 

Theoretical,  numerical,  and  practical  considerations  limit  the  range  of  wave  numbers 
used.  Scales  consistent  with  typical  MM5  use  are  adopted  so  that  Ax  A z.  As  noted  by 
Leutbecher  and  Volkert  (2000),  using  a  wavelength  of  two  horizontal  grid  cells  would  result 
in  large  errors,  due  to  the  artificial  dissipation,  which  is  designed  to  suppress  such  waves. 
Therefore  wavelengths  are  selected  so  that  at  least  four  horizontal  grid  cells  represent  a  wave. 
The  vertical  wavelength  must  be  large  enough  to  be  well  resolved  on  the  given  grid  and  small 
enough  so  that  waves  do  not  span  large  changes  in  background  density  as  required  by  (19) 
for  the  calculation  of  Cg.  For  the  atmosphere  used,  this  is  equivalent  to  an  upper  bound  on 
the  vertical  wavelength  of  about  6000  m.  In  practice,  the  vertical  wavelength  is  kept  well 
below  this  limit  to  maximize  the  time  until  reflections  from  the  upper  and  lower  boundaries 
interfere  with  the  analysis.  For  convenience,  waves  are  also  selected  on  the  basis  of  beam 
angles  and  wave  periods. 

An  example  of  a  beam  that  has  developed  but  has  not  yet  been  influenced  by  boundary 
effects  is  shown  as  a  cross-section  of  the  perturbation  density  field  in  Figure  2.  Such  images 
are  reviewed  at  different  stages  of  beam  evolution  to  confirm  the  flow  is  two-dimensional, 
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Case 

Horizontal  spacing  (m) 

Vertical  spacing  (m) 

^ X 

period  (min) 

Cgz/Cgx 

1 

3000 

434 

.000265 

60 

.11 

2 

3000 

217 

.000228 

60 

.13 

3 

2500 

434 

.00024 

60 

.13 

4 

2700 

217 

.00027 

60 

.12 

5 

3000 

434 

.000154 

120 

.053 

6 

3000 

217 

.00013 

120 

.065 

7 

2500 

434 

.000174 

120 

.072 

8 

2700 

217 

.000134 

120 

.064 

Table  2:  MM5  grid  spacing  and  waves. 


not  impacted  by  boundary  reflections,  and  that  the  waves  satisfy  the  dispersion  relation. 

The  eight  MM5  cases  are  summarized  in  Table  2.  The  cases  represent  waves  with  either  a 
60-minute  period  or  with  a  120-minute  period,  simulated  by  four  different  grids.  The  waves 
are  subject  to  filtering  at  the  boundary  as  they  enter  as  well  as  to  numerical  dissipation  as 
they  propagate.  Therefore,  input  wave  parameters  may  differ  from  observed  parameters  in 
the  MM5  results.  For  calculation  purposes,  two  observable  wave  parameters  are  taken  from 
these  results  to  characterize  the  waves.  The  wave  period  observed  was  consistent  with  its 
input  value,  so  this  was  used  as  one  parameter.  For  the  other  parameter,  the  horizontal  wave 
number,  kx  was  used,  since  it  could  be  measured  in  a  consistent  way  from  case  to  case,  using 
images  like  the  one  in  Figure  2. 

The  vertical  wave  number  is  calculated  using  the  numerical  dispersion  relation.  This  was 
accomplished  using  the  forward  problem  that  solves  for  frequency,  or  period,  and  decay  rate 
(|A|),  given  horizontal  and  vertical  wave  number  and  grid  parameters.  The  procedure  was 
applied  repeatedly  to  obtain  the  inverse,  that  is  to  obtain  the  vertical  wave  number,  given 
the  horizontal  wave  number  and  the  frequency.  The  results  are  shown  in  Table  3,  along  with 
the  decay  rate  and  the  slope  of  the  group  velocity  for  each  of  the  eight  cases  based  on  the 
numerical  analysis. 


5  Results 


The  beams  produced  for  the  eight  cases  are  evaluated  both  in  the  context  of  dissipation 
associated  with  viscous  theory  and  dissipation  expected  based  on  numerical  analysis.  Viscous 
theory  provides  a  gauge  for  the  levels  of  dissipation  observed  in  the  numerical  experiments. 
Numerical  analysis  provides  a  means  to  predict  how  numerical  dissipation  varies  with  wave 
and  grid  parameters. 

For  comparison  to  viscous  theory,  MM5  beam  amplitudes  as  a  function  of  the  beam- 
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Case 

kz 

|A| 

Cgz/C'gx 

1 

-.00120 

.9968 

.24 

2 

-.00106 

.9982 

.22 

3 

-.00110 

.990 

.23 

4 

-.00250 

.9978 

.23 

5 

-.000154 

.9996 

.086 

6 

-.00124 

.9998 

.11 

7 

-.00161 

.9997 

.12 

8 

-.00128 

.9998 

.11 

Table  3:  Calculated  wave  parameters  based  on  numerical  analysis 


following  coordinate  £,  are  fit  to  the  function 

=  (95) 

to  find  the  exponent  F.  The  value  of  F  is  compared  to  the  the  value  2/3  from  (22)  to 
indicate  the  MM5  level  of  damping  relative  to  theoretical,  viscous  damping.  An  example 
of  amplitude  along  a  beam  is  shown  in  Figure  3  for  case  3  (from  Table  3),  with  markers 
to  indicate  MM5  values,  a  solid  curve  to  show  the  best-fit  function  for  those  values,  and 
a  dashed  line  to  show  theory.  In  this  example,  MM5  amplitude  decreased  faster  than  the 
theoretical  prediction  for  a  beam  in  a  viscous  fluid.  The  oscillations  in  the  MM5  values  are 
largely  due  to  the  beam  center  following  a  staircase  pattern  through  the  MM5  grid. 

A  summary  of  F  values  for  all  eight  cases  is  shown  in  Figure  4.  This  shows  that  for  cases 
1-4,  MM5  damping  is  greater  than  theoretical  viscous  damping,  while  for  the  remainder 
it  is  less.  This  suggests  that  in  cases  where  kxAx  >  .55,  or  when  fewer  than  about  eleven 
horizontal  grid  cells  represent  a  wave,  MM5  numerical  dissipation  is  stronger  than  what 
would  be  expected  due  to  molecular  viscosity,  while  in  cases  with  more  than  about  eleven 
horizontal  grid  cells  per  wave,  it  is  weaker.  Numerical  analysis  of  dissipation  is  used  to 
extend  results  over  a  range  of  waves  and  grids. 

In  order  to  show  that  numerical  dissipation  analysis  does  in  fact  describe  MM5  perfor¬ 
mance,  the  two  are  compared.  First,  the  results  of  the  numerical  experiment  are  quantified 
in  terms  of  amplification  factor,  which  leads  to  dissipation  value.  The  experimental  ampli¬ 
fication  factor  is  based  on  fluid  speeds  observed  along  a  beam,  V (£).  Values  of  V (£)  at  a 
fixed  time  are  interpreted  as  equivalent  local  speeds  at  different  times.  For  example,  losses 
over  the  period  At,  when  waves  travel  from  £0  to  £i,  are  given  by  the  amplification  factor 

)"  =  (96) 

where  N  =  (£i  —  Co)/ {CgAt)  is  the  number  of  large  time  steps.  In  this  formulation  the  square 
root  of  density  is  included  because  V ^/po  is  conserved  along  a  compressible,  inviscid  beam. 
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£(kin) 


Figure  3:  Amplitude  along  a  beam  from  MM5  at  a  single  time,  best  fit  curve,  and  theory 
based  on  case  3  in  Table  3. 


Figure  4:  Beam  amplitude  loss  in  terms  of  F  for  the  eight  numerical  experiment  cases,  along 
with  theoretical  value,  F  =  2/3. 
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-4 


Figure  5:  Comparison  of  experimental  dissipation  {Ve  +  Ba)  with  dissipation  predicted  by 
analysis  (Va  +  Ba)  for  the  eight  numerical  experiment  cases. 

The  amplification  factor  due  to  numerical  experiment,  Ae,  is  computed  as  the  least  squares 
fit  value  for  Av  to  all  beam  positions  from  a  single,  upper  beam.  Typically,  ten  to  twenty 
grid  points  along  a  beam  are  available,  depending  on  grid  spacing  and  other  parameters. 

Following  the  form  of  the  dissipation  according  to  analysis  given  by  Equation  (94),  a 
provisional  dissipation  value  based  on  experiment  is  In  (1  —  Ae ).  This  can  be  written  as  the 
sum  of  a  constant  and  a  term  proportional  to  the  order  of  dissipation,  denoted  Ve.  To  unify 
the  constants  added  to  Ve  and  Va ,  the  average  difference  between  the  constants  is  used  as 
an  adjustment  term,  assuming  Da  approximates  Ve.  The  final  form  for  the  experimental 
dissipation  value  is 

Ve  +  Ba  =  In  (1  -  Ae)  -  (in  (1  -  Ae )  -  In  (1  -  Aa))  ,  (97) 

where  the  bars  indicate  averages  over  the  eight  cases  studied. 

The  experimental  dissipation,  Ve+Ba,  is  compared  to  calculations  of  Va+Ba  for  the  eight 
simulated  beams  in  Figure  5,  using  case  numbers  from  Table  2  as  markers.  The  diagonal  line 
indicates  perfect  agreement  between  analysis  and  experiment,  with  the  amount  of  dissipation 
increasing  toward  the  top  right  of  the  diagonal.  Markers  are  close  to  the  line,  indicating 
the  analysis  is  an  approximate  predictor  of  numerical  dissipation.  The  agreement  between 
experiment  and  analysis  shown  in  Figure  5  allows  for  conclusions  that  apply  beyond  the 
experimental  cases,  using  the  numerical  dissipation  analysis  in  a  limited  range  of  parameters. 

The  results  on  Figure  5  form  two  clusters:  an  upper  cluster  of  the  short-wave  cases 
numbered  1-4  and  the  lower  cluster  of  the  long-wave  cases  numbered  5-8.  Of  the  two 
clusters  the  short-wave  cluster  has  greater  numerical  dissipation  as  indicated  by  its  higher 
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Figure  6:  Beam  amplitude  predicted  by  analysis  for  a  range  of  horizontal  and  vertical  pa¬ 
rameters. 


position  on  the  diagonal.  This  is  consistent  with  Figure  4  in  which  the  short-wave  cases 
are  shown  to  contain  greater  dissipation  than  the  long-wave  cases.  Given  that  the  short 
and  long  waves  differ  by  a  factor  of  two,  while  the  horizontal  grids  differ  by  only  about  ten 
percent,  it  is  not  surprising  that  results  are  partitioned  by  wavelength  and  not  horizontal 
grid  spacing.  However,  the  vertical  grid  spacings  differ  by  a  factor  of  two,  but  do  not  produce 
anywhere  near  the  impact  of  the  same  factor  when  applied  to  horizontal  wavelength.  For 
example,  cases  1  and  2,  between  which  vertical  grid  spacing  differs  by  a  factor  of  two,  show 
little  difference  in  dissipation  in  comparison  to  the  dissipation  difference  between  cases  1  and 
5,  which  reflect  a  factor-of-two  difference  in  horizontal  wavelength.  The  same  observation 
can  be  made  with  other  pairings  of  the  data  in  Figure  5.  This  confirms  the  importance  of 
horizontal  grid  spacing  noted  earlier  by  Leutbecher  and  Volkert  (2000). 

Figure  6  shows  the  predicted  amplitude  based  on  the  amplification  factor  determined  by 
the  numerical  analysis  of  Section  3.  Figure  6  shows  that  amplitude  increases,  or  equivalently 
dissipation  decreases,  as  either  kxAx  or  kzAz  decreases.  In  other  words  numerical  dissipation 
decreases  with  refined  resolution  of  waves.  The  dissipation  is  so  much  more  sensitive  to 
horizontal  refinement  than  vertical  refinement  is  evidenced  by  the  shape  of  the  surface  in 
Figure  6.  The  values  (kx,  kz)  =  (0.00036,  —0.0033)  are  used  for  the  figure,  with  uo  satisfying 
the  dispersion  relation  and  At  =  3Ax/1000.  Similar  results  are  found  for  other  wavenumber 
pairs. 
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6  Conclusions 


Idealized  experiments  demonstrate  that  MM5  simulation  of  gravity-wave  beam  propa¬ 
gation  is  consistent  with  theoretical  predictions  of  speed  and  direction.  MM5  beams  lose 
amplitude  due  to  numerical  dissipation.  Using  molecular  viscosity  as  a  gauge,  the  numer¬ 
ical  dissipation  can  be  greater  or  less  than  viscous  dissipation  predicted  for  a  compressible 
atmosphere,  depending  on  the  grid  resolution  and  wave  parameters. 

As  expected,  improving  the  resolution  of  the  wave  in  either  the  horizontal  or  vertical 
generally  decreases  numerical  dissipation.  Both  numerical  analysis  and  numerical  experiment 
show  greater  sensitivity  to  horizontal  resolution  than  vertical,  reflecting  the  influence  of  the 
horizontally-based  artificial  dissipation. 

Due  to  the  relative  sensitivity  to  horizontal  spacing,  as  compared  to  vertical  spacing,  the 
horizontal  spacing  alone  can  be  used  to  approximately  predict  numerical  dissipation.  The 
amount  of  numerical  dissipation  for  a  wave  represented  by  eleven  horizontal  grid  cells  is 
approximately  equivalent  to  the  amount  of  dissipation  that  would  be  imposed  if  molecular 
viscosity  were  present.  Numerical  dissipation  is  greater  than  this  amount  if  the  wave  is 
resolved  by  fewer  horizontal  cells,  and  less  than  this  amount  if  more  horizontal  resolution  is 
used. 

These  results  could  be  used  to  better  understand  MM5  results  and  as  a  consideration  in 
selecting  grids.  This  approach  establishes  a  way  to  perform  idealized  numerical  experiments 
that  can  be  applied  to  other  models  for  model  investigation  or  to  compare  different  models 
on  the  basis  of  their  core  calculations  as  opposed  to  their  parameterizations  or  boundary 
conditions. 

Further  work  is  in  order  to  make  use  of  the  predictions  of  numerical  analysis  in  terms  of 
wave  quantities. 
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